﻿using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
using MathNet.Numerics.Providers.FourierTransform;
using System.Numerics;
using UtilZ.Dotnet.Ex.Model;

namespace UtilZ.Dotnet.Ex.FourierTransform
{
    unsafe public partial class FourierEx
    {
        private static void PrimitiveForwardReal(in SpanZ<float> samples, SpanZ<Complex32> temp, FourierOptions scaling)
        {
            for (int i = 0; i < temp.Length; i++)
            {
                temp[i] = new Complex32(samples[i], 0.0f);
            }

            Forward(temp, scaling);

            samples[0] = temp[0].Real;
            samples[1] = 0f;
            for (int i = 1, j = 2; i < temp.Length / 2; i++)
            {
                samples[j++] = temp[i].Real;
                samples[j++] = temp[i].Imaginary;
            }

            var n = temp.Length;
            if (n.IsEven())
            {
                samples[n] = temp[temp.Length / 2].Real;
                samples[n + 1] = 0f;
            }
            else
            {
                samples[n - 1] = temp[temp.Length / 2].Real;
                samples[n] = temp[temp.Length / 2].Imaginary;
            }
        }

        private static void PrimitiveForwardReal(in SpanZ<double> samples, SpanZ<Complex> temp, FourierOptions scaling)
        {
            for (int i = 0; i < temp.Length; i++)
            {
                temp[i] = new Complex(samples[i], 0.0f);
            }

            Forward(temp, scaling);

            samples[0] = temp[0].Real;
            samples[1] = 0f;
            for (int i = 1, j = 2; i < temp.Length / 2; i++)
            {
                samples[j++] = temp[i].Real;
                samples[j++] = temp[i].Imaginary;
            }

            var n = temp.Length;
            if (n.IsEven())
            {
                samples[n] = temp[temp.Length / 2].Real;
                samples[n + 1] = 0f;
            }
            else
            {
                samples[n - 1] = temp[temp.Length / 2].Real;
                samples[n] = temp[temp.Length / 2].Imaginary;
            }
        }









        private static void BackwardReal(in SpanZ<float> spectrum, SpanZ<Complex32> temp, FourierTransformScaling scaling)
        {
            var n = temp.Length;
            temp[0] = new Complex32(spectrum[0], 0f);
            for (int i = 1, j = 2; i < temp.Length / 2; i++)
            {
                temp[i] = new Complex32(spectrum[j++], spectrum[j++]);
                temp[temp.Length - i] = temp[i].Conjugate();
            }
            if (n.IsEven())
            {
                temp[temp.Length / 2] = new Complex32(spectrum[n], 0f);
            }
            else
            {
                temp[temp.Length / 2] = new Complex32(spectrum[n - 1], spectrum[n]);
                temp[temp.Length / 2 + 1] = temp[temp.Length / 2].Conjugate();
            }

            Backward(temp, scaling);

            for (int i = 0; i < temp.Length; i++)
            {
                spectrum[i] = temp[i].Real;
            }
            spectrum[n] = 0f;
        }

        private static void Backward(SpanZ<Complex32> spectrum, FourierTransformScaling scaling)
        {
            if (spectrum.Length.IsPowerOfTwo())
            {
                if (spectrum.Length >= 1024)
                {
                    Radix2InverseParallel(spectrum);
                }
                else
                {
                    Radix2Inverse(spectrum);
                }
            }
            else
            {
                BluesteinInverse(spectrum);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                    {
                        HalfRescale(spectrum);
                        break;
                    }
                case FourierTransformScaling.BackwardScaling:
                    {
                        FullRescale(spectrum);
                        break;
                    }
            }
        }


        private static void BackwardReal(in SpanZ<double> spectrum, SpanZ<Complex> temp, FourierTransformScaling scaling)
        {
            var n = temp.Length;
            temp[0] = new Complex(spectrum[0], 0f);
            for (int i = 1, j = 2; i < temp.Length / 2; i++)
            {
                temp[i] = new Complex(spectrum[j++], spectrum[j++]);
                temp[temp.Length - i] = temp[i].Conjugate();
            }
            if (n.IsEven())
            {
                temp[temp.Length / 2] = new Complex(spectrum[n], 0f);
            }
            else
            {
                temp[temp.Length / 2] = new Complex(spectrum[n - 1], spectrum[n]);
                temp[temp.Length / 2 + 1] = temp[temp.Length / 2].Conjugate();
            }

            Backward(temp, scaling);

            for (int i = 0; i < temp.Length; i++)
            {
                spectrum[i] = temp[i].Real;
            }
            spectrum[n] = 0f;
        }

        private static void Backward(SpanZ<Complex> spectrum, FourierTransformScaling scaling)
        {
            if (spectrum.Length.IsPowerOfTwo())
            {
                if (spectrum.Length >= 1024)
                {
                    Radix2InverseParallel(spectrum);
                }
                else
                {
                    Radix2Inverse(spectrum);
                }
            }
            else
            {
                BluesteinInverse(spectrum);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                    {
                        HalfRescale(spectrum);
                        break;
                    }
                case FourierTransformScaling.BackwardScaling:
                    {
                        FullRescale(spectrum);
                        break;
                    }
            }
        }
    
    
    }
}
